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ABSTRACT 

We study several versions of the Schmidt-Kennicutt (SK) relation obtained for isolated 
spiral galaxies in TreeSPH simulations run with the GADGET3 code including the novel 
MUlti-Phase Particle Integrator (MUPPl) algorithm for star formation and stellar feedback. 
This is based on a sub-resolution multi-phase treatment of gas particles, where star formation 
is explicitly related to molecular gas, and the fraction of gas in the molecular phase is com- 
puted from hydrodynamical pressure, following a phenomenological correlation. No chemical 
evolution is included in this version of the code. The standard SK relation between surface 
densities of cold (neutral-nmolecular) gas and star formation rate of simulated galaxies shows 
a steepening at low gas surface densities, starting from a knee whose position depends on disc 
gas fraction: for more gas-rich discs the steepening takes place at higher surface densities. Be- 
cause gas fraction and metallicity are typically related, this environmental dependence mimics 
the predictions of models where the formation of H2 is modulated by metallicity. The cold 
gas surface density at which HI and molecular gas surface densities equate can range from 
~ 10 up to 34 M0 pc^^. As expected, the SK relation obtained using molecular gas shows 
much smaller variations among simulations. We find that disc pressure is not well represented 
by the classical external pressure of a disc in vertical hydrostatic equilibrium. Instead is well 
fit by the expression Pfit = ScoidCcoidK/6, where the three quantities on the right-hand side 
are cold gas surface density, vertical velocity dispersion and epicyclic frequency. When the 
"dynamical" SK relation, i.e. the relation that uses gas surface density divided by orbital time, 
is considered, we find that all of our simulations stay on the same relation. We interpret this 
as a manifestation of the equilibrium between energy injection and dissipation in stationary 
galaxy discs, when energetic feedback is effective and pressure is represented by the expres- 
sion given above. These findings further support the idea that a realistic model of the structure 
of galaxy discs should take into account energy injection by SNe. 
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1 INTRODUCTION 

In the last decade a significant step forward in the phenomenologi- 
cal understanding of star formation in galaxies has been achieved, 
thanks to many observational campaigns of nearby (see, e.g., 
IBoissier et al.|2007t Kennicutt et al. 2007] [Walter et al.|2008t and 
distant (e.g. [Bouche et al.||2007[ |Daddi et al.||2010[ [Genzel et aT] 
|2010[ l galaxies. In particular, the relation between surface densi- 
ties of gas and star formation rate (hereafter SFR), the so-called 
Schmidt-Kennicutt relation jSchmidt| 19591 |Kennicutt|1998| here- 
after standard SK), has acquired a higher degree of complexity, 
passing from a simple power law to a much more structured and 



environment-dependent relation. Here by environment we mean the 
local properties, averaged over ~k;pc scale, of a galaxy patch that 
hosts star-forming molecular clouds. We will call standard, molec- 
ular, HI and dynamical SK the relations obtained by putting on 
the y-axis the surface density of SFR Esfr, and on the x-axis, re- 
spectively, the total cold gas surface density Scold, the molecu- 
lar gas surface density Emoi, the HI gas surface density Ej// or 
the total cold gas surface density divided by the orbital time-scale, 
Ecoid/Torb, where Torb ~ 2Txr/V{r) (V(r) being the galaxy rota- 
tion curve). 

Wide consensus has recently been reached on the idea that the 
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standard SK is a reflection of the more "fundamental" molecular 
SK ( |Wong & Blitz|2002| [Blitz & Rosolow sky 200 4i[2006) . while 
the HI SK is weak if not absent jBigieTeral.|2008|[2010^ . Indeed, 
in normal (and non-barred) spiral galaxies the fraction of gas in 
molecular clouds increases toward the galaxy center, while the HI 
gas surface density, T,hi, saturates at a value of ~ 10 pc^^ 
and declines in the inner kpc. The old notion of a star formation 
threshold in disc galaxies (Martin & Kenn icutt|2001| l has thus been 
revised to a steepening of the standard SK at low surface densities 
( [Boissier et al.[[2007{ [Bigiel et al.[[2008] l, driven by the declining 
molecular fraction. 

The gas surface density at which the transition from HI- to 
molecular-dominated gas takes place, or equivalently the satura- 
tion value of T,Hi, has been proposed to be a function of galactic 
environment ([Krumholz et al.|2009b||Gnedin & Kra vtsov[20l"0l[Pa^ 
[padopoulos & P elupessy 2010, see also Schaye 2004J, with dwarf 
galaxies like the Magellanic Clouds showing higher values (Bo- 
latto et al. 2009; [Bolatto et al.|201I| see also references in [Fuma-j 
[galli et al.|2010| ). This is in line with the high-redshift evidence of 
low efficiency of star formation in Damped Lyman-alpha systems 
at 2 ~ 2 ( [Wolfe & Chen|20 06 1, that are thought to be the external, 
gas-rich, metal-poor regions of young disc galaxies. 

The slope of the molecular SK relation is debated. For 
THINGS spiral galaxies, [Bigiel et al.[ ( [2008[ l report a slope of 
1.0 ± 0.2 when measured on a spatial grid of 750 pc bin size. An 
average value of ~1 has been confirmed very recently by [Bigiel] 
[et al.[ ^201 \) . A steeper slope, more consistent with the canonical 
1.4 value of [Kennicutt|1998[ is reported by [Kennicutt et al.[p007) 
for star-forming regions in M5 la. [Llu et al.| ( [20 1 1 1 interpret this dis- 
crepancy as an effect of subtraction of background emission in Ho? 
and dust, and claim that a super-linear slope results when proper 
subtraction is performed. At higher surface densities, a steeper re- 
lation is suggested by observations of Ultra Luminous Infra-Red 
Galaxies (ULIRGs) and Sub-mm Galaxies (SMGs) 
2007| l, but recent observations (Da ddi et ar][2010[ 
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20T0| [Boquien et aL][20Tl ) suggest that at 2 ~ 2 ULIRGS and 



SMGs, on the one side, and non-IR bright and BzK galaxies, on 
the other side, trace two parallel molecular SK relations with slope 
~ 1.4 and separated by nearly one order of magnitude in normal- 
ization (IR-bright galaxies having higher Esfr). 

The interpretation of this phenomenological picture is still 
under discussion. Based on observations, the declining molecu- 
lar fraction with gas surface density was proposed by [Blitz &[ 
[Rosolowsky[ p004[ [2006| ) to be driven by external pressure, i.e. 
the midplane pressure of gas in vertical hydrostatic equilibrium in 
a thin disk. However, this relation has large scatter and is as scat- 
tered as other relations with, e.g., disc mass surface density (jL eroyj 
et al.|[2008|. Alternativ ely, many authors (Pelupessy et al. 2006, 
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proposed models where the molecular fraction is regulated by the 
equilibrium balance between production of H2 and destruction of 
the same molecule by UV radiation from young stars (an assump- 
tion recently criticized by |Mac Low & Glover|20I0) . Both creation 
and destruction channels are regulated by dust abundance, because 
dust is both a catalyst and a shield. As a consequence, the molec- 
ular fraction is predicted to be driven by gas surface density (or 
column density) and modulated by gas metallicity. Therefore, the 
scaling of molecular fraction with gas surface density, or equiva- 
lently the saturation value of Tjhi, should be a function of metal- 
licity. |Fumag^iret^|2010| have recently tested the two assump- 
tions (pressure-driven or gas surface density-driven molecular frac- 



tion) against data on nearby spiral and dwarf galaxies, and report a 
marginal preference for the second hypothesis. 

The varying slope of the of the molecular SK, steepening to- 
ward high Ecoid from ~ 1.0 to ~ 1.4 — 1.7, has been interpreted 
by [Krumholz et al.[p009bl l as an effect of the decoupling of molec- 
ular clouds in normal spiral galaxies from the rest of the Inter- 
Stellar Medium (ISM). These authors argue that molecular clouds 
are known to have a roughly constant surface density and pressure 
(and then dynamical time, see , Solomon et al.[1987p , so they are not 
in pressure equilibrium with the rest of the ISM and their consump- 
tion time Emoi/Esfr results to be ~ 2 Gyr ( Bigiel et al.|2008 1, ir- 
respective of the molecular gas surface density computed on ~kpc 
scale. This last quantity is indeed to be considered as a measure 
of the filling factor of molecular clouds. The decoupling breaks at 
higher gas surface densities, where the ISM is able to pressurize 
the molecular clouds so that their dynamical time scales again with 
the inverse of the square root of the density at ~kpc scales. In this 
regime the molecular SK takes values more similar to the canonical 
1.4 one. 

The evidence of a double standard (or molecular) SK at high 
redshift has only been proposed very recently. [Teyssier et al.[j20l"0l > 
have presented a hydrodynamic simulation, performed with the 
AMR RAMSES code ( [Teyssier|2002^ , of two merging galaxies re- 
sembling the Antennae. They found that, when the force resolu- 
tion reaches values as small as 12pc, the predicted standard SK is 
boosted and reaches the relation found by Dad di et al.|(2010| l to 
hold for ULIRGs and SMGs. However, the authors do not show a 
simulation, run at the same resolution, of the equivalent of a BzK 
galaxy that lies on a relation one order of magnitude lower. 

A second way of expressing a "star formation law" is by cor- 
relating the surface density of SFR with Ecom divided by the gas 
orbital time torh- This dynamical SK was suggested by jWyse &[ 
|Siik| ( [T989) , [Slik| ( [T997l > and [Elmegreen| ( [T997l > to be more "funda- 
mental" than the standard SK, on the basis of the influence that disc 
rotation and shearing exert on star-forming clouds. It has attracted 
less attention than the standard SK, and in most cases it has been 
reported as equally acceptable from the observational point of view 
(e.g. [Kennicuttj 1 998] > . More recently, [Tan[ ^2010^ tested against ob- 
servations of local galaxies the hypothesis of a linear standard SK 
compared to several other proposals of star formation "laws". He 
found that data do not favour this linear hipothesis. At higher red- 
shift, Daddi et al.[ l |201 0T) noticed that the dicothomy in the molec- 
ular SK of normal and IR-bright galaxies disappears when the dy- 
namical SK is considered. 

The nature of the SK relation and its environment-dependent 
nature is of fundamental importance in the field of galaxy forma- 
tion, because SK-like relations are widely used in models to predict 
the amount of stars formed in a star-forming galaxy. In particular, in 
many N-body hydrodynamic simulations (see, e.g. [Katz et al.|I996l > 
the SFR of gas particles or cells is computed as e x Mgas/rdyn, 
where e is an efficiency and the dynamical time is computed on the 
average density of the particle. In many of the most recent simla- 
tions of galaxy formation it is assumed that the SFR obeys by con- 
struction a standard SK law with a cut at low density (e.g. [Springel[ 
[& Hemquist1[2003l l. [ScFaye & Dalla Vecchiaj ( |205g) showed that 
this is equivalent to assuming an effective equation of state for star- 
forming particles of the kind P oc p^"^^ . Only higher-resolution 
simulations, resolving scales below 50 pc, are able to follow the 
formation of mol ecules ^Pelupessy et alT][2 006, Pelupessy & Pa^ 
[padopoulos[2009[[Robertson & Kravtsov|2008)|Gnedin et al.|2009l l 
and then to predict molecular fractions and kpc-averaged SK re- 
lations, but the cost of these simulations makes it very hard with 
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presently available facilities to push even one of them to z = in 
a full cosmological context. 

One exception in this sense is given by the MUPPI model 
(MUlti-Phase Particle Integrator; [Murante et al.||2010l hereafter 
paper I), based on a sub-resolution multi-phase treatment of the 
star-forming gas particles and developed within the GADGET 
TreePM-l-SPH code (Springel 20051. In this model, the star forma- 
tion rate is not forced to follow a SK-like relation, so its adherence 
to this law is a prediction of the model. The details of the model 
are described below, but two points are worth mentioning. First, in- 
spired by the observational result of [Blitz & Rosolowsky| ( (2006| l, 
the molecular fraction of the cold component of a gas particle is 
scaled with the hydrodynamic pressure of the SPH particle. Sec- 
ond, the multi-phase treatment allows thermal feedback from SNe 
to efficiently heat the gas. In paper I, free parameters were tuned to 
reproduce the observed SK relation in the case of an isolated spiral. 
Milky Way-like galaxy. We noticed that the SK relation traced by 
an isolated low surface brightness dwarf spiral galaxy differs from 
the Milky Way one in the same way as spirals and dwarfs differ in 
the data of [Bigiel et aLjpOOSt . 

In this paper we show the SK relations resulting from a set of 
MUPPI simulations of isolated halos, including the ones used in pa- 
per I. The value of this analysis is not only to present the results of a 
specific model but also to understand how the SK relations depend 
on galactic environment when the disc is heated by feedback and 
the molecular fraction is scaled with pressure and does not depend 
on metallicity. Section 2 describes the MUPPI model and the initial 
conditions used for our simulations. Section 3 presents the result- 
ing SK relations. Interpretation of results requires an assessment of 
the vertical structure of simulated discs, presented in Section 3.2. 
Section 3.3 is devoted to a discussion of the dynamical SK relation. 
Section 4 gives a discussion of present results in comparison with 
available literature, while Section 5 presents the conclusions. 



2 SIMULATIONS 

The main properties of the simulated galaxies used in this paper are 
reported in Table 1 . 

The first set of initial conditions is the one used in paper I, 
and were generated following the procedure described in [Springel] 
( |2005| l. They are near-equilibrium distributions of particles consist- 
ing of a rotationally supported disc of gas and stars and a dark 
matter halo. For MW and MW_HR a stellar bulge component is 
included. Bulge and halo components are modeled as spheres with 
[Hemquist] ( |1990| l profiles, while the gaseous and stellar discs are 
modeled with exponential surface density profiles. To start from a 
relaxed and stable configuration, we first evolve the galaxy mod- 
els for 10 dynamical times with non-radiative hydrodynamics. We 
then use the configurations evolved after 4 dynamical times as ini- 
tial conditions for our simulations. To test the effect of resolution 
on the SK relations, the MW galaxy is used also at a higher resolu- 
tion (MW_HR). 

The second set of initial conditions has been used by Spri ngel] 
[& Hemquist[ l [2003[ l for a resolution test of their star formation 
code. Gas particles ar e embedded in a 1.39 x 10^° Mg static NFW 
I Navarro et aL]]l996 l| halo (lO"' h''^ Mg with h = 0.72), and 
rotate at a speed corresponding to a spin parameter of A = 0.1, 
radially distributed following |Bullock et al.[ ( |2001| l. Gas is initially 
in virial equilibrium with the halo. When cooling is switched on, 
gas particles slowly coalesce into a rotating disc. With this sim- 
ple setting it is possible to use very high force resolution, so that 



the vertical structure of the disc is well resolved. These initial con- 
ditions are available at four resolutions. We show results for only 
the highest one. We checked that results are very stable when the 
resolution is degraded. 

In all cases we forced the SPH smoothing length of gas parti- 
cles not to drop below 1/2 of the Plummer-equivalent softening. 

2.1 The model 

MUPPI has been developed within a non-public version of the 
GADGET2 Tree-PM-l-SPH code (Springel 2005) that includes an 
entropy-conserving formulation of SPH (^Springel & Hemquist[ 
[2002| l. It has already been ported into the more efficient GADGET3 
code. All details are given in [Murante et al.[ ( |2010[ >, while we only 
describe here some relevant features of the MUPPI algorithm. 

This model is inspired by the multi-phase analytic model of 
star formation and feedback by jMonaco [ ( |2004b^ . A particle enters 
the multi-phase regime when its temperature is lower than a thresh- 
old (set to 5 X 10* K) and its density, recast in terms of particle 
number density (with a molecular weight of ^ = 0.6) is higher than 
0.01 cm~'^. This threshold is an order of magnitude lower than the 
commonly used value of ~ 0.1 cm"'', which is typically tuned to 
obtain a cut in the standard SK relation at a gas density ~ 10 Mg 
pc-^ 

A multi-phase particle is assumed to be made up of three com- 
ponents: two gas phases and a stellar phase. The two gas phases are 
assumed to be in thermal pressure equilibrium. As in [Springel &[ 
[Hemquist]j2003^ , the cold gas phase is assumed to have a tempera- 
ture Tc = 1000 K, while the temperature of the hot gas phase is set 
by the particle entropy. Upon entrance into the multi-phase regime, 
all the mass is assumed to reside in the hot phase; cooling does not 
lead to a lowering of the temperature but to a deposition of mass in 
the cold phase. A fraction of the cold mass is assumed to be in the 
molecular form. [Blitz & Rosolowsky| (2006 l (see also jLeroy et al.[ 
]2008j ) found that the ratio between molecular and HI gas surface 
densities correlates with the so-called external pressure, which is 
the pressure expected at the galaxy midplane in the case of a thin 
disc composed by gas and stars in vertical hydrostatic equilibrium. 
This was estimated using a simplified version of the expression pro- 
posed by [Elmegreen[p989y : 



fcxt — -^GEcold (Scold + 



(1) 



Here R = a^oid / c"* is the ration between the vertical r.m.s. velocity 
dispersions of cold gas and stars (here a denotes velocity dispersion 
while E denotes surface density). The exponent a of the correlation 
Emoi/Sm oc P^t was found to be 0.9 ± 0.1. Inspired by this 
finding, and adopting an exponent of 1 for simplicity, we use the 
following equation to estimate the molecular fraction: 



/mol(P) 



1 + Po/P 



(2) 



Here P is the hydrodynamical pressure of the SPH particle (dif- 
ferent from the external pressure used in the observati onal cor- 
relation), and we adopt Po/k = 35000 K cm"'' as in 
Rosolowsky] ( (200 6)l^ 



Blitz & 



The three components exchange mass through four mass 
flows: cooling, star formation, restoration and evaporation. Cooling 

^ We have checked that our simulations produce relations similar to the ob- 
servational one (though with significant scatter) when pressure is estimated 
as in [Blitz & Rosolowsky| l [2006} . 
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Figure 1. Gas surface density (left column: face-on; middle column: edge-on) and temperature (right column) maps of simulated galaxies. Color coding is 
reported in the color bars, numbers refer to the Log of surface density or temperature. The galaxies shown are, from top to bottom, MW, MW_HR, DW and 
SH. 
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Table 1. Basic characteristics of simulated galaxies. Column 1: Simulation name. Coluirm 2: Gravitational Plummer-equivalent (P-e) softening for gas 
particles. Column 3: Mass of DM halo. Column 4: Mass of DM particle. Column 5: Stellar mass. Column 6: Mass of star particle. Column 7: Half-mass radius 
of stars. Column 8: Cold gas mass. Column 9: Initial mass of gas particle (before spawning stars). Column 10: Half-mass radius of cold gas. Column 11: Gas 
fraction. Notes (1): for MW, MW_HR and DW stellar masses include both the disc and bulge (only MW and MW_HR) stars present in the initial conditions 
and the newly formed stars, which are a minority. (2): in the above cases we report the stellar mass particle of old stars, because new stars (whose particle mass 
is TTlgas /4) give a negligible contribution to the disc. (3): The DM halo in the SH simulation is static. 
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(Mq) 


(Mq) 


(kpc) 


fraction 


MW 


0.69 


9.4 ■ 10" 


3.5 ■ 10*^ 


4.2 ■ 10"' 


1.3 ■ 10^ 


4.8 


3.3 ■ 10^ 


7.4 ■ 10-* 


5.6 


7.3% 


MW_HR 


0.41 


9.4 • 10" 


6.9 • 10'^ 


4.2 ■ 10"' 


2.6 ■ 10^ 


4.4 


3.2 • 10^ 


1.5 ■ lO-i 


5.4 


7.1% 


DW 


0.42 


1.6 • 10" 


8.1 • 10^ 


7.8 ■ 10^ 


1.6 ■ 10^ 


8.5 


1.9 • 10^ 


3.9 ■ 10-* 


8.3 


20% 


SH 


0.042 


1.4 ■ 10"' 


_(3) 


1.4 ■ 10^ 


2.2 ■ 10^ 


0.77 


1.4 ■ 10^ 


8.7- 10^ 


5.2 


99% 



is computed with a standard cooling function assuming zero metal- 
licity. Thermal energy resides in the hot phase, so cooling is com- 
puted using its density, which is much lower than the average one. 
Star formation takes place in the molecular phase, with a consump- 
tion time-scale proportional to the particle dynamical time tdyn: 

M* = /mol(P)/*A/cold/tdy„(nc) (3) 

Here /moi is the pressure-dependent molecular fraction of equa- 
tion |2] /* = 0.02 is a parameter of star formation efficiency, de- 
termining the fraction of a molecular cloud that is transformed into 
stars per dynamical time, and Afcoid the mass of the cold phase 
in the particle. The dynamical time tdyn is computed on the cold 
phase as soon as this collects 90 per cent of the total gas mass, 
and is frozen until the particle exits from the multi-phase regime 
(see paper I for a detailed discussion of this hypothesis). A thing 
is worth noticing in equation |3] because of the hypothesis of pres- 
sure equilibrium and constant temperature of the cold phase, ric is 
proportional to pressure P; at the same time, the fraction of gas 
mass in the hot phase is always very low, so /moi is very similar to 
the fraction of total gas in molecular form. As a consequence, the 
particle star formation rate is primarily regulated by gas pressure 
(with the complication that tdyn is computed at the beginning of a 
star formation cycle and then kept frozen, while /moi is computed 
at each time-step). 

The SFR term of equationjsjdeposits a fraction (1 — frc) of the 
transformed mass into the stellar component, while a fraction /rc, 
restored from massive stars in an Instantaneous Recycling Approx- 
imation (IRA), is given back to the hot phase. The formed stellar 
component is accumulated within the particle, and contributes to its 
inertia but not to its gas mass in all SPH computations. The evap- 
oration rate is assumed to be due to the destruction of molecular 
clouds and amounts to a fraction /cv = 0.1 of the SFR. 

Production of star particles is done according to the stochas- 
tic star formation algorithm of |Springel & Hemquist| ( [2003| > (see 
paper I for details). We allow for 4 generations of star particles to 
be spawned by each parent gas particle. Each new star particle is 
produced at the expense of the stellar component and, if needed, of 
the cold phase. 

The three components of a multi-phase particle are of course 
assumed to be subject to the same hydrodynamical forces. To alle- 
viate the effect of this unphysical assumption, and to mimic the 
destruction of a star-forming cloud after a few dynamical times 
( |Monaco|2004a^ , the code forces each particle to leave the multi- 
phase regime after two dynamical times tdyn, computed as speci- 
fied above. 

One SN is generated each Af*,sjv = 120 Mq of stars formed. 



and each SN generates lO^'^ erg of energy. Of the energy generated 
in the IRA, a small fraction /fb,i = 0.02 is given to the local hot 
phase to sustain its high temperature, while a fraction /fb,o ~ 0.3 
is distributed to the hot phases of neighbour particles in a 60-degree 
wide cone anti-aligned with the gas density gradient. Energy con- 
tributions to particles are weighted by their distance from the cone 
axis, to mimic the expansion of SN-driven blasts along the least re- 
sistance path ( |McKee & Ostriker|I977) . The present version of the 
code distributes only thermal energy. 

The assumption of a molecular fraction regulated by pressure 
(equation |2]l is very important because it makes the evolution of 
the system intrinsically runaway: star formation generates SNe, en- 
ergy feedback from SNe pressurizes the hot phase, the increase in 
pressure leads to an increase in molecular fraction and thus to an 
increase in SFR. The runaway halts when the molecular fraction 
saturates to unity. However, the dynamical response of the pressur- 
ized particle is able to limit this runaway through the expansion 
work done on neighbours. This intrinsic runaway behaviour, to- 
gether with the long cooling times, are the main reasons for our 
efficient thermal feedback. 



3 THE SK RELATIONS IN SIMULATED DISCS 

The first three simulations, the MW (at standard and high resolu- 
tion) and DW test cases of paper I, start from already formed disc 
galaxies with 10 and 20 per cent gas fractions, and are shown at 
0.5 Gyr, after the initial transients due to the switcliing on of stel- 
lar feedback have (almost) died out and while gas consumption is 
still negligible. The fourth simulation, the SH halo, forms in an 
isolated, static halo filled with rotating gas initially in virial equi- 
librium; we analyze the simulation at 0.7 Gyr, i.e. at the peak of its 
SFR, when the disc is still largely gas-dominated. In all cases, con- 
clusions are unchanged when simulations are considered at other 
times. We show in figure[T]face-on and edge-on maps of gas surface 
densities and temperatures for the four simulations. It is interesting 
to notice that in the regions interested by star formation the discs 
(most of the MW and MW_HR discs and the inner few kpc of SH) 
are relatively hot and surrounded by thick coronae of gas heated 
by feedback and circulating above or below the disc in a galactic 
fountain. The effect of star formation on the DW galaxy is much 
less evident. 

3.1 The standard, HI and molecular SK relations 

Figure |2] shows the standard, HI and molecular SK relations of the 
four simulations, compared with the data of [Bigiel et al.| p008| > 
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Figure 2. SK relations for simulated galaxies. The black contour levels re- 
port the data from Bigiel et al. (1998), binned in bins of size 0.5 dex; lev- 
els correspond to 1, 2, 5 and 10 observational points per bin, gas surface 
densities include helium. Triansles and thick lines sive results for the four 
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Figure 3. Molecular fraction as a function of cold gas surface density for 
the four simulations. The dotted line inarks the 1/2 value, where HI and 
inolecular gas surface densities are equal. 



for normal spiral galaxies. Gas surface densities are always meant 
to include contribution from helium. In the upper panel the thin 
line represents the fit proposed by |Kennicut"t| jl998^ . Simulations 
have been processed as follows. Our analyses are restricted to cold 
gas; in this paper by cold gas we mean the cold phase of multi- 
phase particles plus all single-phase particles colder than 10 ' K. 
The molecular gas surface density is computed using the molecular 
fraction of equation|2] applied only to multi-phase particles; HI gas 
is just cold minus molecular gas. A galaxy frame is defined by the 
inertia tensor of stars and cold gas, the z-axis corresponding to the 
largest eigenvalue. The angular momentum of the same particles is 
always found to be at most a few degrees off the z-axis. Then, radial 
surface densities (in cylindrical coordinates) of (cold, HI, molecu- 
lar) gas and SFR are computed; these are reported in the figure as 
colored thick lines. The same quantities are computed on a square 
grid in the x-y plane, with bin size of 750 pc, as in the |Bigiel et al.| 
(2008 1 paper; these are reported as triangles, with the same color as 
the corresponding thick lines. 

The following points are apparent from the figure: (i) the sim- 
ulations trace different SK relations when the total or HI gas are 
used, (ii) These different relations correspond to different transi- 
tions from HI- to molecular-dominated gas, or equivalently to dif- 
ferent saturation values ofTiHi- (iii) The differences go in the di- 
rection highlighte d by|Bigiel et al.|p008) , [Bolatto| ( |2009t ; |Bolatto| 
|et al.| \2Q\l) and [Bigiel et al.| ( |2010| l of a higher saturation Y.hi 
in dwarf galaxies, (iv) The molecular SK relation is much tighter; 
also, it is steeper than the Bigiel et al. ( 2008 1 one and consistent 
with |Liu et al.| ( |20l"T| >. (v) The scatter around the SK relation is gen- 
erally smaller than the data, especially for the SH simulation, (vi) 
As for the MW and MW_HR simulations, the SK relation is rather 
stable with resolution; this is confirmed by applying the same anal- 
ysis to the SH halo at lower resolution. 

To better quantify the difference between the simulations, in 
figure |3] we show the relation between surface density and molec- 
ular fraction /moi, a quantity that is directly comparable with data. 
The Ecoid value at which the molecular fraction is 1 /2 ranges from 
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Figure 5. Correlation of the ratio of tlie radial average of Psim/fcxt with 
Toomre parameter Qtot for the four simulations. Color coding is reported 
in the panel; thick continuous lines correspond to that part of the galaxy 
beyond two softening lengths from the center and with significant SFR, 
thin dashed lines refer to the rest of the galaxy. The dashed thick line has 
slope 1/3. 

10 M0 pc"^ of MW to 35 M0 pc^^ of SH. This is in agreement 
with the measures of the Large Magellanic Cloud, though more 
extreme cases like the Small Magellanic Cloud are not recovered 
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As pointed out by |Fumagalli et al.]j2010^ , if the molecular SK 
is the "fundamental" relation then a pressure law for the molecular 
fraction results in an environment-dependent standard SK. This is 
true because, as evident in equation [T| for pressure in the case of 
vertical hydrostatic equilibrium, this quantity depends not on gas 
surface density alone (the quantity in the x-axis of the standard SK) 
but on gas and stellar surface densities, the latter being multiplied 
by the ratio R of gas and star vertical velocity dispersions (see also 
|Shietal.|20lT) . In other words, the cut in the standard SK is deter- 
mined by the (velocity-weighted) gas fraction. So our result is easy 
to interpret as long as hydrodynamic pressure of our discs scales 
with gas fraction in a similar way as the external one. It is worth 
noticing at this point that our simulated galaxies have gas fractions 
ranging from low to very high values, so they span most of the in- 
teresting range for this quantity. 

However, the validity of equation [T| obtained in the case of 
gas in vertical hydrostatical equilibrium, must be checked in our 
simulations where energy is continually injected in discs. This is 
done in next sub-section. 

3.2 Vertical structure of simulated discs 

In figure [4] we report, for the four simulated galaxies, the relation 
between pressure and Ecow as found in the simulations. We show 
in each panel the hydrodynamical pressure found in simulations, 
Psim, averaged in bins of the x-y grid as colored triangles, and the 
radial profile of the average of the same quantity as a line with 
a darker color. The external pressure Pcxt is computed using equa- 



Figure 6. Gas vertical velocity dispersions, including thermal and kinetic 
energy J4j, for the four simulations. Color coding is reported in the panel. 

tion^from the radial profiles (colored dot-dashed lines). Moreover, 
to ease the comparison, each panel reports, as grey thick lines, the 
radial Ecoid — Psim relations from the other simulations. To obtain 
Pcxt, vertical velocity dispersions ) of cold gas and star parti- 
cles have been computed in the galaxy frame. For the stars, this 
quantity is equated to aj, while for the gas we use: 

l^coM = (I'z) + , (4) 

where Cs is the gas sound speed, computed on the average parti- 
cle temperature and density. Figure [4] shows that the relation be- 
tween Psim and Ecoid varies in a very similar way as that between 
Pcxt and Ecoid- This confirms the validity of our interpretation for 
the environmental dependence of the standard SK. However, exter- 
nal and simulated pressures show significant dicrepancies that are 
worth addressing. 

These discrepancies may be related to the fact that we are 
comparing the expected midplane pressure with the average one 
over the whole disc height. While for barely resolved discs the two 
quantities cannot differ much, the SH simulation allows us to test 
to what extent midplane and average pressures are comparable. As 
a matter of fact, particles in the midplane show a broad range of 
pressures, because the accumulation of cold mass at the beginning 
of a star formation cycle causes a depressurization by one order of 
magnitude, while successive star formation re-pressurizes the star- 
forming particles. As a result, pressure does not show a smooth 
decrease with height on the disc, and the (mass-weighted) average 
of pressure over the disc height is always very similar to the mid- 
plane one. A similar trend was found by |Tasker & Bryan] pOOSl l, 
using a much better resolution, when energy feedback from SNe is 
considered. 

We noticed that in our simulations the ratio of true pressure 
Psim and Pcxt, computed the first in radial bins and the second 
from radial profiles, correlates well with the Toomre parameter Q 
of the disc, computed at the same radius. For a disc made of a single 
component with surface density E and velocity dispersion a, we 
have Q(r) = ctk/ttGE, where k = l/(r)^2 + 2d\nV/dlnr/r) 
is the epicyclic frequency and V{r) the rotation curve. Being our 
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Figure 4. Relation between total cold gas surface density and hydrodynamical pressure in simulations. Triangles denote the relation in 750pc bins, while darker 
continuous lines of the same color give the radial averages. The bright dot-dashed and dashed lines give radial estimates based respectively on hydrostatic 
equilibrium (equation[TJ and on equation|6] Each panel contains the results from one simulation, while only the radial averages of hydrodynamical pressure 
from the other simulations are reported in each panel as grey lines. 



discs composed by stars and gas, we compute the disc total Toomre 
parameter Qtot by using the simple approximation of |Wang & Silk| 
|1994|, that is correct for our discs characterized by relatively high 
velocity dispersions l |Romeo & Wiegert|201H : 



1^1 

?cold Q* 



fcoldK 



7rG(Ecold + 



(5) 



It is easy to show that Qtot = o-coidKEcoid/2Poxt- 

Figure |5] shows for our four simulations and at all radii the re- 
lation between Qtot and the ratio Psim/-Pext- The most interesting 
regions with E^fr > 10~^ M© yr~^ kpc~^ and radius larger than 
two softening lengths have been highlighted as thick continuous 



lines. The two quantities show a linear correlation that is well fit by 
Psim/fcxt = Qtot/3; the round number 3 adapts well to the SH 
simulation where the vertical structure of the disc is fully resolved. 
Then, a much better fit to the Psim — Scow relation is given by: 



Pflt = Poxt X 



6 



.IdCcoldK • 



(6) 



The second equality is obtained using equations [T] and [5] Pflt is 
shown in figure|4]as dashed lines; it gives a much better approxima- 
tion to the Psim — Ecoid relation, with some discrepancy for the DW 
galaxy at E^oid < 5 M0 pc~^, where Esa is anyway very low. We 
checked the validity of Pflt by comparing it to the pressure profile 
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Figure 7. Toomre Q parameters for cold gas (dotted lines), stars (dashed 
lines) and total (continuous lines). Color coding is reported in the panel. 



of many versions of our galaxies, obtained in our tests by varying 
model parameters and assumptions. In particular, we tested changes 
in the computation of dynamical time used in equation|3]in several 
ways, e.g. by equating it to the average one at entrance into multi- 
phase or by recomputing it at each time-step. We found in all cases 
that equation[6]always gives a good fit to pressure, with significant 
discrepancies found only at very low pressure. Therefore we con- 
sider this result to be independent of details of our sub-resolution 
model. Presently, we have no analytical explanation of why multi- 
plying the external pressure by Qtot /3 gives a better fit to the av- 
erage SPH pressure in our simulations. We interpret equation |6] as 
the quasi-equilibrium pressure of a disc with continuous injection 
of energy from SNe: hotter discs, with high Qtot, are characterized 
by higher pressure than Pcxt, while marginally stable discs have a 
lower pressure by a factor up to ~ 3 (a similar result was found by 
[Koyama & Ostriker|2009[ see section 4 for more details). 

More insight on the vertical structure of the disc can be ob- 
tained by analysing the quantities (TcoM and Qtot, that enter in the 
computation of Pflt . 

Figure |6] shows the values of gas vertical velocity dispersions 
fcoid for the four galaxies, computed both on the 750 pc grid and 
in radial profiles. Clearly the assumption of a constant velocity dis- 
persion, often made in the literature to represent this quantity, is a 
poor approximation of our results. Moreover, these velocities are 
much higher than the ~ 5 — 10 km s^^ usually assumed for galaxy 
discs[3 

The values of the Q parameters for the four simulations are 
shown in figure |7] Qtot was computed as in equation [5] the 
same parameter was computed only for gas and stars as Qi — 
aiK/nCEi (where i is either cold or ★), and its values are re- 
ported in the figure with dotted and dashed lines, respectively. The 
MW and DW discs have Qtot ~ 1 in the central, star-forming re- 



^ This quantity includes contribution from the thermal energy of the hot 
phase, so it is not directly comparable to observations. We will return on 
this in section 4. 



gions, while Q, assumes higher values. Qcow is very high, ^ 20; 
this quantity cannot be directly compared with observations of cold 
gas, because the main contribution to (Tcoid comes from the sound 
speed, that is determined by the thermal energy of the hot phase 
of multi-phase particles. If the sound speed is neglected, we obtain 
Qcoid ~ 5 — 10. With this caveat, these values of the Toomre pa- 
rameters are in very good agreement with what found by |Leroy| 
|et al.| poos ') in nearby galaxies, where Qtot is just above unity 
while Q*, and especially Qcoid take up higher values. At the same 
time, the gas-dominated SH disc is much more stable and kinemat- 
ically hot. As a conclusion, while our stellar-dominated discs tend 
to regulate to keep Qtot ~ 1, this is not a general rule. 

As a further test, we run our simulations, starting from the 
same four sets of initial conditions but using the effective model of 
|Springel & Hemquist| ( |2003) , that is known not to provide efficient 
thermal feedback. We obtained colder discs, with higher Qtot pa- 
rameters and hydrodynamical pressure not well fit either by Pcxt or 
by our equation|6] The resulting standard SK relation is in line with 
that originally fo und by |Springel & Hemquist|p003| l: it lies on the 



Kennicutt 



jl998| relation above ~5 Mq pc , and cuts below that 



threshold, not following the mild steepening found in data. This 
surface density threshold somehow depends on the galaxy: we find 
it at 3, 4.5 and 6 M© pc"^ for the MW, DW and SH simulations. 
This dependence is easy to understand: the cut is determined by the 
imposed volume density threshold for star formation. Since disc 
temperatures do not vary much (especially in these colder discs) 
the relation between volume and surface density of gas scales sim- 
ilarly to that between pressure and surface density. Of course no 
prediction can be made in this case for the HI and molecular SK. 
We conclude that a simple volume density threshold for star for- 
mation, used in conjunction with an ineffective feedback scheme, 
gives a standard SK with environmental dependence that goes in 
the same direction as the one presented in this paper, but cannot 
reproduce the data at the same level of detail. 



3.3 The dynamical SK 

Figure [8] shows the dynamical SK of our four simulations. Unfor- 
tunately, a dataset like that of ,Bigiel et al.| ( (2008) is unavailable for 
this relation; most observations give global estimates of galaxies. 
As a consequence, we do not compare this relation with data in this 
paper. 

Remarkably, the four simulations, that produced different 
standard SK's, now trace a unique, non-linear relation, with a scat- 
ter that is even lower than that in the molecular SK. To make this 
equality more evident, we decrease the point size of the grid-based 
relation to let radial profiles be more visible. 

It is easy to show, using equations [T] and [5] that a "universal" 
dynamical SK relation, valid for all galaxies, can be obtained under 
the following simple assumptions: (i) the disc is in vertical hydro- 
static equilibrium, so that pressure is well represented by Pcxt; (ii) 
the disc is marginally stable, with a Qtot Toomre parameter equal 
to 1; (iii) the velocity dispersion of gas acoid is constant for all 
galaxies; (vi) the rotation curve is flat so that k ~ v^/Torb; (v) 
SFR is a function of pressure. But we have shown above that the 
first three conditions don't hold for our discs, so this simple inter- 
pretation cannot be valid here. On the other hand, we noticed that 
our simulated discs are rather hot when we take into account the 
effective temperature of gas particles, which is determined by the 
thermal energy of the hot phase. In the following we demonstrate 
that the dynamical SK can be interpreted as the result of the bal- 
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Figure 8. Dynamical SK relations for simulated galaxies. Points and thick 
lines give results for the four simulations in 750pc bins and in radial profiles 
in cylindrical coordinates; color coding is given in each panel. Point size has 
been decreased to make the agreement of the four relations more evident. 



ance between energy injection and dissipation, without assuming 
that SFR is directly influenced by disc rotation. 

Energy is continually injected by SN feedback and lost by 
radiative cooling and viscous dissipation. The equilibrium among 
these processes can be illustrated with the help of a simplified ap- 
proach that quantifies the energy made available to gas particles, 
in both kinetic and thermal form, after cooling has radiated away a 
part of it. Injection of kinetic and thermal energy can be expressed 
as (see also equation 14 of paper I): 



(7) 



where the constant v^^ = (/fb,i + /fb,o)10^^ erg/Af*.s„ takes ac- 
count of the energy made available for feedback by MUPPI and e 
is an efficiency parameter that quantifies the fraction of injected 
energy that is not radiated away by cooling. This injected energy 
is then transformed by expansion into kinetic energy and is later 
dissipated by viscosit}]^ this keeps the disc in a quasi-stationary 
state. For a disc of height H^s, disc energy will likely be dissipated 
on one disc height crossing time, icrosa = ^^off/o"coid- Notably, 
this is the same rate at which turbulence is dissipated, as long as 
the driving length of turbulence is the size of the typical SN-driven 
bubble (e.g. |Mac Low|2003] l, which must be ~ J/gff if bubbles die 
by blowing out of the disc (see also| Monaco|2004bt . 

Let's define the surface density of energy in the disc as: 



ri _ 3„ 2 

■ijE — -Zjcoldfcold J 



(8) 



where we assume equipartition between the three translational de- 
grees of freedom. Then the energy dissipation rate is: 



^ The fact that in our simulations viscosity is the numerical one included in 
our SPH code is at this stage immaterial, as long as this behaves similarly 
to dissipation of turbulence that would take place if the resolution were 
adequate to describe it. 
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In the definition of effective disc height HbR = Scoid/2pcoid 
the midplane density can be substituted with pressure using the 
equation of state P — PcoidC^oid- Using equation [6] for the pres- 
sure, it is easy to show that: 
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It then follows that: 
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In a stationary system, the effective energy injection will 
equate dissipation, so that Einj = Ediap- This implies that: 



or, using the approximated value of «: ~ \/2/ 
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This relation allows us to recast the interpretation of a unique dy- 
namical SK relation in terms of energy injection: it will hold as 
long as, at fixed Estr, the post-cooling efficiency of energy injec- 
tion e scales with the square of the gas vertical velocity dispersion, 
i.e. the gas specific energy. This can be seen in the other direction: 
the specific energy of the gas f^^oid rritist scale with the fraction of 
specific thermal energy that has time to perform PdV work before 
cooling radiates it away, ev'i^^. This is a property that arises naturally 
from our sub-resolution feedback model, at least for values of free 
parameters that do not widely differ from the fiducial ones selected 
in paper I; as in the case of pressure, we tested the validity of this 
result with a large suite of MUPPI simulations on the same sets of 
initial conditions, and many combinations of parameters and phys- 
ical assumptions. The result of an environment-dependent standard 
SK and a unique dynamical SK holds in all cases where a disc effi- 
ciently heated by feedback and in a stationary state is obtained. 

According to our interpretation, the fact that our four simula- 
tions all lie on the same dynamical SK is no evidence that SFR is 
directly determined by galaxy rotation, at variance with the moti- 
vation that has been used to introduce this relation. To make this 
more clear, the dispersion rate of energy given in equation|9]can be 
written, without loss of generality, as Edisp = SPcTcoid- Pressure 
in our simulations is well reproduced by equation |6] that depends 
on the kinematical state of the disc through Qtot, and then on k. 
As a result, the dynamical stationarity of the disc induces a depen- 
dence of SFR on the epicyclic frequency, i.e. on the orbital time. 
So, the fact that our four simulations all lie on the same dynamical 
SK is telling us that they are stationary discs kept out of vertical 
hydrostatic equilibrium by SN feedback, and not that their SFR is 
directly determined by rotation or shearing - it is only indirectly in- 
fluenced by it through the dependence of pressure on Qtot- Indeed 
we noticed that some of our simulated discs stay out of this dynam- 
ical SK, and this happens when, e.g., a bar instability takes place. In 
this case the condition of stationarity is violated until a new equil- 
brium configuration is reached. So these discrepancies confirm our 
interpretation that simulations trace a unique dynamical SK as long 
as they are in a stationary condition. 

When run with the effective model of |Springel & Hernquist| 
([2003), the three galaxies (MW, DW, SH) show similar but not 
imique dynamical SK, the differences being comparable to those 
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in the standard SK (section [3^ . Tliis is a result of inefficient tlier- 
mal feedback in this model: as long as pressure is not well fit by 
equation [6] and the injection of energy is marginal, we do not ex- 
pect galaxies to lie on a unique dynamical SK. 



4 DISCUSSION 

One result of this paper is that our pressure-based model of star for- 
mation produces a standard SK with a knee at a gas surface density 
that depends on gas fraction in a way that resembles the metallicity 
dependence proposed by |Krumholz et al.| ( |2009b) and [Gnedin &| 
IKravtsov] ^20 1 0^ . If the molecular fraction is regulated by pressure, 
gas-rich galaxies become dominated by molecular gas at higher gas 
surface densities. It is easy to show that the gas surface density at 
which Smoi equates Effj, which we call Eeq, scales with the gas 
fraction fi = Ecoid/(Scoid + E*), as: 

E ^ (14) 

Here Pq is the normalization of equation |2] In most chemical evo- 
lution models, /i is related to metallicity; for instance, in the simple 
closed-box model, Z — Y ln(l//i), where Y is the metal yield per 
stellar generation. Then, a pressure-based molecular fraction can 
mimic a metallicity dependence. But the maximum value of Ecq 
is reached for n = 1, in which case Eoq.max = y^2Po/TrG ~ 
34 Mq pc~^; this is the Eeq value of the SH simulation, that has a 
99 per cent gas fraction. According to |Bolatto| ( |2009^ , Eeq for the 
Large Magellanic Cloud is similar to 34 Mq pc~ , while for the 
Small Magellanic Cloud it is ~100 Mq pc~^. This is confirmed 
by a quick comparison with Krumholtz et al's model: an increase 
in Ecq by a factor of 3, as large as the difference between MW 
and SH, is obtained by a decrease of metallicity by a similar factor, 
which is relatively modest. The same conclusion can be reached by 
considering the saturation Ehi value, that reaches ~ 20 Mq pc~^ 
for our SH disc while, for instance, [Fumagalli et al.| ( [20I0[ l report 
higher values for a few low metallicity dwarf galaxies. We con- 
clude that a pressure-driven molecular fraction cannot explain the 
whole observed range of variation of Eeq. Because the motivation 
for molecular fraction being directly modulated by metallicity is 
very strong (but see the comments by |Mac Low & Glover|20I0l l, it 
is well conceivable to construct mixed scenarios where equation |2] 
is valid and the normalization Pq depends on metallicity. 

The search for a physical interpretation of our results lead us to 
an expression for gas pressure in star-forming discs given by equa- 
tion [6] Discs subject to this pressure have some interesting prop- 
erties: for instance, neither E* nor a* appear in the expression of 
pressure (though the gravity of stars enters in determining acow), 
gas effective height is related to the ratio of gas velocity dispersion 
and epicyclic frequency (equation 1 10^ and the sound crossing time 
of effective height is simply proportional to the orbital time (equa- 
tion[TTJ. Comparison with data is not straightforward, as a^oM in- 
cludes a major contribution from the sound speed of the hot phase 
(we give further comments below) and direct pressure estimates are 
hard to obtain. Nonetheless, these predictions can in principle be 
tested against observations of the Milky Way and nearby galaxies, 
so they may constitute a basis for a theory of the non-equilibrium 
vertical structure of discs effectively heated by feedback. 

Our simulated discs show total vertical velocity dispersions 
(figure[6} that are well in excess of the ~ 10 km s~^ value that is 
usually assumed to hold for discs. However, CTcoW in our simula- 
tions is dominated by the thermal sound speed (equation |4j com- 



puted at the particle effective temperature, and this last quantity is 
determined by the thermal content of the hot phase. This means 
that (Jcoid cannot be compared with the velocity dispersion of cold 
clouds in real galaxies. At the same time, vertical velocity disper- 
sions of gas particles (computed without the sound speed) of MW 
and DW simulations show values that are in much better agreement 
with those measured in THINGS galaxies by |Tamburro et al.|p009| ) 
(see figures 1 1 and 18 in paper I). In our simulations, multi-phase 
particles are composed by two phases at the sub-grid level, but are 
seen as a single entity by the SPH code. This means that, as long as 
a multi-phase cycle goes on, the hot phase is not free to leave the 
disc, while the cold phase is pulled by the former. This artificial en- 
trainment results, at the macroscopic level, in a velocity dispersion 
of gas particles that is realistic and, from a detailed comparison of 
MW and MW_HR, rather stable with resolution in that part of the 
disc that is well resolved in both simulations. This means that our 
effective model is correctly producing a gas disc that is thermally 
warm but kinematically colder. It must be noticed that the hot coro- 
nae we produce around the discs have temperatures ~ 3 x lO*" K 
and densities ~ 10~^ cm"'', so their presence is likely ruled out 
by X-ray observations (e.g. |Crain et aLpOl O*). However, insertion 
of metal cooling would change this energy balance in favour of 
kinetic energy, so we expect this hot corona to be less prominent 
when chemical evolution is properly taken into account. 

Our results thus suggest that a complete modeling of a disc 
heated by feedback must fully take into account the multi-phase 
nature of the ISM, where the ~kpc-height corona of warm/hot gas 
surrounding a spiral galaxy may have an important dynamical role. 
A step forward in the modeling of a spiral disc subject to con- 
tiuous production and dissipation of energy was recently taken by 
|Elmegreenl ( |20I l| l, who addressed the stability of a disc where en- 
ergy is continually injected and dissipated over some multiple of 
the crossing time of turbulence. However, in that paper only radial 
and tangential perturbations were considered and vertical equilib- 
rium was assumed. 

The results presented in this paper rely on how well the ver- 
tical structure of the disc is resolved in simulations. In the MW 
and DW cases effective disc heights are of the same order of the 
gravitational softening (that was kept of the same order as the one 
used in cosmological simulations with the same mass resolution), 
so numerical convergence of the results should be demonstrated. 
We showed results for the MWJTR simulation, and found no clear 
dependence on resolution in any of our results. Moreover, the light 
disc forming out of the SH simulation has an effective height of 
~ 1 kpc, with a gravitational softening of only 43 pc, so in this 
case the vertical structure is well resolved. We conclude that, de- 
spite the vertical structure of the discs is barely resolved in some 
of our simulations, the results should be not strongly affected by 
resolution. 

Due to the difficulty in performing measures of gas pressure, 
observers typically use equation[T]to estimate pressure (often mak- 
ing strong assumptions on velocity dispersions), but its validity has 
rarely been tested. [Koyama & Ostriker ( 2009 1 compared this for- 
mula with simulations of turbulent ISM in a shearing disc. In their 
simulations the gas disc is assumed to be much thinner than the 
stellar one, heating terms take account of cosmic rays. X-rays and 
H2 formation and destruction, while radiative feedback from mas- 
sive stars is modeled as localized increases of heating rate, but no 
SNe are present. Their spatial resolution is of order of ~1 pc. They 
found that midplane and mass-weighted pressures typically differ 
by an order of magnitude, and that equation [T| is very close to the 
armonic mean of the two. Interestingly, their midplane pressure is a 



12 P. Monaco et al. 



factor of three lower than Poxt, which is what we find for Qtot = 1- 
Our results are not comparable to their simulation, due to the vastly 
different resolution and to the much hotter state of our discs (they 
have (Jcoid ~ 5 km s^^). They also proposed an improved analytic 
estimate of midplane pressure; we compared it with our Psi-m and 
found that it does not improve much with respect to Pext • Our sim- 
ulations are instead comparable to those of |Tasker & Bryan] j2008^ 
and Joung et al. ( 2009 ), that have resolutions of 25 and 2 pc respec- 
tively, include feedback frm SNe and show multi-phase structures 
of the ISM in broad agreement with the assumptions made to de- 
sign MUPPI. Unfortunately, they don't explicitly address the ques- 
tion whether their pressure is well reproduced by the external one 
of |Elmegreeri| ( [T989) . 



5 CONCLUSIONS 

In this paper we have shown how simulations based on the MUlti- 
Phase Particle Integrator (MUPPI) model developed within the 
GADGETS TreePM-l-SPH code, are able to give predictions on the 
various SK relations discussed in the literature. MUPPI is based on 
the assumption, suggested by observations of Blitz & Rosolows^ 
( |2004 2006), that the molecular fraction is modulated by pressure. 
Moreover, it has the interesting feature of making thermal feedback 
effective and able to heat discs. So, the interest of this paper does 
not rely only in the test of a specific sub-resoution model for star 
formation but it allows to understand what are the testable conse- 
quences of a pressure-based molecular fraction in discs efficiently 
heated by SN feedback. 

Our main conclusions are the following. 

(i) A pressure-based molecular fraction produces an 
environmental-dependent standard SK relation, owing to the fact 
that the relation between pressure and gas surface density is mod- 
ulated by gas fraction. This variation is very similar to that found 
between spiral and dwarf galaxies ("Bigiel et al.|2008||2010^ , and it 
could be interpreted as a metallicity dependence, since gas fraction 
is typically related to metallicity in most chemical evolution sce- 
narios. However, the variation in the quantity Eeq, the gas density 
at which the molecular fraction is 1/2, cannot be larger than a factor 
of 3 between normal spirals and gas-dominated discs, and values of 
Eeq > 34 M© pc^^ cannot be obtained in this framework. 

(ii) We analyzed in detail the vertical structure of our sim- 
ulated galaxies, and found that hydrodynamical pressure is not 
well recovered by the vertical hydrostatic equilibrium value of 
|Elmegreeii| ( |198 9'), because kinematically hotter discs show higher 
pressure. A better fit is given by: 

-Pflt ~ TrEcoldCcoldK 

5 

(equation |6j, that we interpret as the pressure of a disc with con- 
tinuous energy injection. This expression allows to connect the ef- 
fective disc height with the gas velocity dispersion (computed in- 
cluding kinetic and thermal energies) as Hcs ~ ScTcoid/'^ (equa- 
tion[l0j. 

(iii) Quite interestingly, our four simulated galaxies lie on the 
same (non-linear) dynamical SK relation, independently of their 
gas fraction. This is not a straightforward consequence of our as- 
sumptions, and was not expected. It is worth recalling that a similar 
phenomenology was found by |Daddi et aL] ( |20To) in the very dif- 
ferent context of 2 ~ 2 star-forming galaxies. We interpret this 
result as a manifestation of balance between energy injection from 



SNe and energy dissipation. We have shown that, under the hypoth- 
esis that gas energy is dissipated both by cooling and by viscosity, 
and that the latter works on the timescale of one sound crossing 
time of the disc effective height, we can obtain for stationary discs 
a unique dynamical SK if the efficiency of energy injection after 
cooling scales with gas specific energy. This is found to result from 
energy balance in our multi-phase particles under a wide range of 
cases. 

(iv) Other results are of some interest. The model follows well 
the standard, molecular and HI SK relations, with a tight molecular 
SK that is a straightforward consequence of the model assumptions. 
This molecular SK has a slope of ~ 1.4, marginally steeper than 
the 1.2 value found by Bi giel et al.| pi308 ) but in agreement with 
|Liu et al.]j201 1^ . The scatter in the simulated relations is small, and 
this may hint that most scatter is due to observational errors, if not 
to putting together galaxies that lie on slightly different relations. 
However, our sub-resolution model gives by construction the aver- 
age properties of the ISM on scales that are similar to the 750 pc 
scale used to bin the data, and this may be the reason for the low 
scatter. 

These simulations show that energy injection by SNe is fun- 
damental in determining the structure of star-forming discs, and 
that the warm/hot phases created by stellar feedback may have an 
important role in disc dynamics. Future observations will need to 
address the issue of directly determining gas pressure in order to 
test whether the usually assumed formula of'Elmegreen (1989) or 
some different estimates, like that provided by equation|6] apply. 
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